Binder 0

Bienvenue dans ce chapitre interactif utilisant ThebeLab !

Exemple de code interactif

Run code

import matplotlib.pyplot as plt
import numpy as np

def plot_sine(frequency=1):
    x = np.linspace(0, 30, 1000)
    y = np.sin(frequency * x)
    plt.figure(figsize=(8, 4))
    plt.plot(x, y)
    plt.title(f'Sine Wave with Frequency {frequency}')
    plt.xlabel('x')
    plt.ylabel('sin(x)')
    plt.show()

plot_sine()

try:
    import dolfinx
    print("dolfinx importé avec succès")
except ImportError as e:
    print(f"Erreur lors de l'importation de dolfinx : {e}")
    import sys
    print(f"Python path : {sys.path}")

dolfinx importé avec succès

Problème torsion d’arbres

Librairies

Run code

# Importation des librairies
import dolfinx  # Module principal de FEniCSx pour le calcul par éléments finis
from dolfinx import mesh, fem, plot, io, default_scalar_type  # Sous-modules FEniCSx essentiels
from dolfinx.fem.petsc import LinearProblem  # Résolution de systèmes linéaires
from mpi4py import MPI  # Interface pour le calcul parallèle
import ufl  # Unified Form Language pour les formulations variationnelles
import numpy as np  # Calcul numérique efficace
import matplotlib.pyplot as plt
import math  # Module de fonctions mathématiques standard de Python
import pyvista as pv  # Visualisation 3D scientifique
import gmsh  # Générateur de maillage 3D
import meshio  # Lecture/écriture de différents formats de maillage
from dolfinx.io import XDMFFile  # Gestion des fichiers XDMF pour données volumineuses
import dolfinx.fem as fem  # Module FEniCSx pour la méthode des éléments finis
import dolfinx.mesh as mesh  # Module FEniCSx pour la gestion avancée des maillages
from petsc4py import PETSc  # Suite de solveurs numériques pour problèmes à grande échelle
import panel as pn
from myst_nb import glue
import ipywidgets as widgets
from IPython.display import display
# Utiliser Panel pour créer un rendu interactif
pn.extension('vtk')

Géométrie

Définition de la géométrie et du maillage

# Initialisation du plotteur PyVista en dehors de la fonction
p = pv.Plotter(notebook=True)
p.set_background("grey")

# Fonction pour créer et visualiser la géométrie
def create_ellipse(rho_x, rho_y, h=1.0, lc=0.02):
    global domain
    gmsh.initialize()
    gmsh.option.setNumber("General.Terminal", 0)
    gmsh.model.add("ellipse")

    # Création de l'ellipse
    ellipse = gmsh.model.occ.addEllipse(0, 0, 0, rho_x, rho_y)
    curve_loop = gmsh.model.occ.addCurveLoop([ellipse])
    surface = gmsh.model.occ.addPlaneSurface([curve_loop])

    # Extrusion pour obtenir un cylindre elliptique
    gmsh.model.occ.extrude([(2, surface)], 0, 0, h)
    gmsh.model.occ.synchronize()

    # Configuration du maillage
    gmsh.option.setNumber("Mesh.CharacteristicLengthMin", lc)
    gmsh.option.setNumber("Mesh.CharacteristicLengthMax", lc)
    gmsh.model.mesh.generate(3)

    # Sauvegarde et conversion du maillage
    gmsh.write("ellipse.msh")
    msh = meshio.read("ellipse.msh")
    meshio.write("ellipse.xdmf", meshio.Mesh(points=msh.points, cells={"tetra": msh.cells_dict.get("tetra", [])}))
    gmsh.finalize()

    # Lecture du maillage avec FEniCSx
    with XDMFFile(MPI.COMM_WORLD, "ellipse.xdmf", "r") as xdmf_file:
        domain = xdmf_file.read_mesh(name="Grid")
        domain.topology.create_connectivity(domain.topology.dim-1, domain.topology.dim)
    
    # Conversion du maillage pour la visualisation avec PyVista
    u_topology, u_cell_types, u_geometry = plot.vtk_mesh(domain)
    u_grid = pv.UnstructuredGrid(u_topology, u_cell_types, u_geometry)

    # Effacer la scène précédente et ajouter le nouveau maillage
    p.clear()
    # Ajout du maillage à la scène de visualisation
    p.add_mesh(u_grid, 
               show_edges=True,  # Affiche les arêtes du maillage
               scalar_bar_args={
                   "title": "u",  # Titre de la barre de couleur
                   "title_font_size": 24,
                   "label_font_size": 22,
                   "shadow": True,
                   "italic": True,
                   "font_family": "arial",
                   "vertical": False  # Orientation horizontale de la barre de couleur
               })

    # Ajout d'un titre à la visualisation
    p.add_text("Cylindre Mesh", font_size=24, color="black", position="upper_edge")

    # Ajout des limites de la boîte englobante
    p.show_bounds(color="black")

    # Ajout des axes de coordonnées
    p.add_axes(color="black")

    # Définition de la couleur de fond
    p.set_background("grey")
    
    # Mise à jour de la scène
    p.show()
    
    return rho_x, rho_y, h
    
# Widgets interactifs pour a, b, et h
rho_x_slider = widgets.FloatSlider(value=0.2, min=0.1, max=1.0, step=0.01, description='Demi-grand axe')
rho_y_slider = widgets.FloatSlider(value=0.2, min=0.1, max=1.0, step=0.01, description='Demi-petit axe')
h_slider = widgets.FloatSlider(value=1.0, min=0.5, max=3.0, step=0.1, description='Hauteur')

# Interface utilisateur et fonction interactive
ui = widgets.VBox([rho_x_slider, rho_y_slider, h_slider])
out = widgets.interactive_output(create_ellipse, {'rho_x': rho_x_slider, 'rho_y': rho_y_slider, 'h': h_slider})

display(ui, out)

# Fonction pour récupérer les valeurs actuelles des curseurs
def get_slider_values():
    current_rho_x = rho_x_slider.value
    current_rho_y = rho_y_slider.value
    current_h = h_slider.value
    return current_rho_x, current_rho_y, current_h

# Exemple d'utilisation : récupération des valeurs
rho_x, rho_y, h = get_slider_values()
rho = np.array([rho_x, rho_y])  # Demi-axes dans les directions x et y

Unhandled exception
Traceback (most recent call last):
  File "/opt/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/aiohttp/web_protocol.py", line 531, in start
    resp, reset = await task
                  ^^^^^^^^^^
  File "/opt/anaconda3/envs/fenicsx-env/lib/python3.12/asyncio/futures.py", line 287, in __await__
    yield self  # This tells Task to wait for completion.
    ^^^^^^^^^^
  File "/opt/anaconda3/envs/fenicsx-env/lib/python3.12/asyncio/tasks.py", line 385, in __wakeup
    future.result()
  File "/opt/anaconda3/envs/fenicsx-env/lib/python3.12/asyncio/futures.py", line 203, in result
    raise self._exception.with_traceback(self._exception_tb)
  File "/opt/anaconda3/envs/fenicsx-env/lib/python3.12/asyncio/tasks.py", line 314, in __step_run_and_handle_result
    result = coro.send(None)
             ^^^^^^^^^^^^^^^
  File "/opt/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/aiohttp/web_protocol.py", line 448, in _handle_request
    assert self._request_handler is not None
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
AssertionError

interactive_panel = pn.pane.VTK(p.ren_win, width=755, height=400)
from myst_nb import glue
glue("img", interactive_panel)

Configuration du problème

Définition de l’espace fonctionnel

# Définition de l'espace fonctionnel pour le problème
V = fem.functionspace(domain, ("Lagrange", 1, (domain.geometry.dim, )))
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)

Définition des frontière du domaine

# Définition des fonctions pour identifier les différentes faces du cylindre
def down_face(x):
    """Identifie la face inférieure du cylindre (z = 0)"""
    return np.isclose(x[2], 0)

def top_face(x):
    """Identifie la face supérieure du cylindre (z = h)"""
    return np.isclose(x[2], h)

def lateral_face(x):
    """
    Identifie la surface latérale du cylindre elliptique
    Utilise l'équation de l'ellipse : x^2/a^2 + y^2/b^2 = 1
    """
    tolerance = 1e-5  # Tolérance pour la comparaison numérique
    return (np.isclose((x[0]**2 / rho[0]**2 + x[1]**2 / rho[1]**2), 1.0, atol=tolerance)) & (0 <= x[2]) & (x[2] <= h)

# Dimension des facettes (2D pour un domaine 3D)
fdim = domain.topology.dim - 1

# Localisation des entités (facettes) correspondant à chaque face
Sigma_l = mesh.locate_entities(domain, fdim, lateral_face)
Sigma_0 = mesh.locate_entities(domain, fdim, down_face)
Sigma_h = mesh.locate_entities(domain, fdim, top_face)

# Préparation des marqueurs de facettes pour les conditions aux limites
marked_facets = np.hstack([Sigma_l, Sigma_0, Sigma_h])
marked_values = np.hstack([
    np.full_like(Sigma_l, 1),  # Marqueur 1 pour la face latérale
    np.full_like(Sigma_0, 2),  # Marqueur 2 pour la face inférieure
    np.full_like(Sigma_h, 3)   # Marqueur 3 pour la face supérieure
])

# Tri et création des tags de maillage
sorted_indices = np.argsort(marked_facets)
facet_tag = dolfinx.mesh.meshtags(domain, fdim, 
                                  marked_facets[sorted_indices], 
                                  marked_values[sorted_indices])

# Définition de la mesure pour les facettes marquées
ds = ufl.Measure('ds', domain=domain, subdomain_data=facet_tag)


# Création d'un dictionnaire pour mapper les noms aux valeurs numériques
surface_ids = {'Sl': 1, 'S0': 2, 'Sh': 3}

Visualisation des frontière du domaine

# Définition d'une fonction pour appliquer des marqueurs aux facettes
def apply_marker(boundary_facets, marker_array, marker_value):
    """
    Applique un marqueur spécifique aux facettes de la frontière.
    
    Args:
    boundary_facets: Indices des facettes de la frontière
    marker_array: Tableau des marqueurs à mettre à jour
    marker_value: Valeur du marqueur à appliquer
    """
    # Filtrer les facettes pour ne garder que celles appartenant au domaine local
    boundary_facets = boundary_facets[boundary_facets < num_cells_local]
    # Appliquer le marqueur
    marker_array[boundary_facets] = marker_value

# Obtenir le nombre de cellules locales (important pour le calcul parallèle)
num_cells_local = domain.topology.index_map(fdim).size_local

# Initialiser les tableaux de marqueurs pour chaque type de frontière
# On crée 4 tableaux : 3 pour chaque face et 1 pour la combinaison
markers = [np.zeros(num_cells_local, dtype=np.int32) for _ in range(4)]

# Appliquer les marqueurs spécifiques à chaque face
apply_marker(Sigma_0, markers[0], 1)  # Face inférieure (z = 0)
apply_marker(Sigma_h, markers[1], 2)  # Face supérieure (z = h)
apply_marker(Sigma_l, markers[2], 3)  # Surface latérale

# Calculer le marqueur combiné
markers[3] = markers[0] + markers[1] + markers[2]

# Ajouter un marqueur spécifique (4) pour les cellules non marquées
# Cela peut être utile pour identifier le volume intérieur, par exemple
markers[3][markers[3] == 0] = 4

# Créer la connectivité topologique pour les facettes
# Ceci est nécessaire pour certaines opérations sur le maillage
domain.topology.create_connectivity(fdim, fdim)

# Obtenir les données de maillage au format compatible avec PyVista
topology, cell_types, x = dolfinx.plot.vtk_mesh(domain, fdim, np.arange(num_cells_local, dtype=np.int32))

# Création du maillage PyVista à partir des données FEniCSx
grid = pv.UnstructuredGrid(topology, cell_types, x)

def add_plot(ax, marker, color, title, threshold_min):
    """
    Fonction pour ajouter des maillages à une sous-fenêtre de visualisation.
    
    Args:
    ax: Axe de la sous-fenêtre PyVista
    marker: Tableau des marqueurs pour les cellules
    color: Couleur pour les cellules filtrées
    title: Titre de la sous-fenêtre
    threshold_min: Valeur minimale pour le filtrage des cellules
    """
    # Mise à jour des données de cellule du maillage avec les marqueurs
    grid.cell_data["Marker"] = marker
    grid.set_active_scalars("Marker")
    
    # Ajout du maillage complet avec les marqueurs
    ax.add_mesh(grid, show_edges=True, color="cyan", scalar_bar_args={
        "title": "Boundary Marker",
        "title_font_size": 24,
        "label_font_size": 22,
        "shadow": True,
        "italic": True,
        "font_family": "arial",
        "vertical": False
    })
    
    # Application d'un filtre basé sur le seuil pour isoler certaines cellules
    grid_filter = grid.threshold(threshold_min, scalars='Marker')
    ax.add_mesh(grid_filter, color=color, show_edges=True)
    
    # Ajout du titre et des axes à la sous-fenêtre
    ax.add_text(title, font_size=12, color="black", position="upper_edge")
    ax.add_axes(color="black")

# Configuration de la visualisation PyVista avec 4 sous-fenêtres (2x2)
pl = pv.Plotter(shape=(2, 2))

# Sous-fenêtre 1 : Affichage de toutes les frontières
pl.subplot(0, 0)
add_plot(pl, markers[3], "red", "All Boundaries", threshold_min=0.5)

# Sous-fenêtre 2 : Affichage de la frontière inférieure
pl.subplot(0, 1)
add_plot(pl, markers[0], "red", "Down Boundary", threshold_min=0.5)

# Sous-fenêtre 3 : Affichage de la frontière supérieure
pl.subplot(1, 0)
add_plot(pl, markers[1], "red", "Top Boundary", threshold_min=1.5)

# Sous-fenêtre 4 : Affichage de la frontière latérale
pl.subplot(1, 1)
add_plot(pl, markers[2], "red", "Lateral Boundary", threshold_min=2.5)

# Configuration finale et affichage
pl.set_background("grey")
pl.show()

# Utiliser Panel pour créer un rendu interactif
pn.extension('vtk')
interactive_panel = pn.pane.VTK(pl.ren_win, width=600, height=400)
from myst_nb import glue
glue("img", interactive_panel)

Conditions aux limites

# Définition des conditions aux limites de Dirichlet (déplacement nul)
u_D = np.array([0, 0, 0], dtype=default_scalar_type)

# Application des conditions aux limites sur chaque face
bc0_S0 = fem.dirichletbc(u_D, fem.locate_dofs_topological(V, fdim, Sigma_0), V)  # Face inférieure
bc0_Sh = fem.dirichletbc(u_D, fem.locate_dofs_topological(V, fdim, Sigma_h), V)  # Face supérieure
bc0_Sl = fem.dirichletbc(u_D, fem.locate_dofs_topological(V, fdim, Sigma_l), V)  # Surface latérale

Propriétés matériau

import ipywidgets as widgets
from IPython.display import display, clear_output
from ipywidgets import HBox, VBox

# Widgets pour les constantes élastiques en GPa
mu_input = widgets.FloatText(value=80, description="μ = ", step=1)
mu_label = widgets.Label(value="GPa")

lambda_input = widgets.FloatText(value=120, description="λ = ", step=1)
lambda_label = widgets.Label(value="GPa")

# Organisation avec les labels d'unité
mu_box = HBox([mu_input, mu_label])
lambda_box = HBox([lambda_input, lambda_label])

# Fonction pour obtenir les constantes élastiques en Pa
def display_elastic_constants(change):
    mu = mu_input.value * 1e9  # Conversion en Pa
    lambda_ = lambda_input.value * 1e9  # Conversion en Pa
    
    # Afficher les valeurs en Pa
    clear_output(wait=True)  # Clear previous output
    
    return mu, lambda_  # Retourner les valeurs en Pa

# Liaison de la fonction d'affichage aux événements de changement de valeur des widgets
mu_input.observe(display_elastic_constants, names='value')
lambda_input.observe(display_elastic_constants, names='value')

# Affichage des widgets
ui_elastic = VBox([mu_box, lambda_box])
display(ui_elastic)

mu, lambda_ = display_elastic_constants('value')
print(f"Module de cisaillement (μ) en Pa: {mu:.2f}")
print(f"Premier paramètre de Lamé (λ) en Pa: {lambda_:.2f}")

Module de cisaillement (μ) en Pa: 80000000000.00
Premier paramètre de Lamé (λ) en Pa: 120000000000.00

Définition des chargements

import ipywidgets as widgets
from IPython.display import display
import numpy as np

def create_angle_widget():
    # Création du widget curseur
    angle_slider = widgets.FloatSlider(
        value=2.0,  # Valeur initiale à 2 degrés
        min=0,
        max=360,
        step=1,
        description='Angle (°):',
        disabled=False,
        continuous_update=True,
        orientation='horizontal',
        readout=True,
        readout_format='.0f',
    )

    # Widget de sortie pour afficher les valeurs
    output = widgets.Output()

    # Fonction pour mettre à jour l'affichage
    def update_angle(change):
        ang_degrees = change['new']
        ang_radians = np.radians(ang_degrees)
        
        with output:
            output.clear_output(wait=True)
            print(f"Angle : {ang_degrees:.0f}° ({ang_radians:.2f} radians)")
        
        return ang_radians

    # Attacher la fonction de mise à jour au widget
    angle_slider.observe(update_angle, names='value')

    # Afficher le widget et la sortie
    display(widgets.VBox([angle_slider, output]))

    # Initialiser l'affichage
    update_angle({'new': angle_slider.value})

    return angle_slider

# Créer et afficher le widget
angle_widget = create_angle_widget()

# Fonction pour obtenir la valeur actuelle en radians
def get_angle_radians():
    return np.radians(angle_widget.value)

# Exemple d'utilisation

# Obtention de l'angle de torsion en radians à partir du widget
ang_radians = get_angle_radians()

# Calcul de alpha (taux de torsion)
# alpha représente le taux de rotation par unité de longueur
# Il est calculé en divisant l'angle de torsion par la longueur du cylindre (h)
# Unités : radians / mètre (rad/m)
alpha = ang_radians / h
print(f"Taux de torsion (α) = {alpha:.4f} rad/m")

# Calcul du rayon équivalent du cylindre elliptique
# On utilise la moyenne quadratique des demi-axes pour obtenir un rayon équivalent
# Cette formule donne une bonne approximation pour les calculs de torsion
R = np.sqrt((rho[0]**2 + rho[1]**2) / 2)
print(f"Rayon équivalent (R) = {R:.4f} m")

# Calcul de la rigidité en torsion C
pi = math.pi  # Nombre Pi
# La formule pour C est basée sur la théorie de la torsion pour un cylindre circulaire
# mu : module de cisaillement du matériau
# R^4 : le rayon à la puissance 4 représente l'influence de la géométrie
# alpha : taux de torsion calculé précédemment
# Le facteur 1/2 est spécifique à la formule de rigidité en torsion
C = mu * pi * (R ** 4) / 2 * alpha
print(f"Rigidité en torsion (C) = {C:.4e} N·m²")

Taux de torsion (α) = 0.0349 rad/m
Rayon équivalent (R) = 0.2000 m
Rigidité en torsion (C) = 7.0184e+06 N·m²

# Calcul de la constante A pour le champ de torsion
A = 2 * C / (pi * R**4)  # A est un scalaire constant positif
print(f"Densité de contrainte de torsion (A_τ) = {A:.4e} N/m³")

# Définition des coordonnées spatiales
x = ufl.SpatialCoordinate(domain)

# Définition des composantes du champ de torsion (vecteur des Contraintes de Cauchy)
T1 = -A * x[1]  # Composante x du champ de torsion / Contrainte de cisaillement dans la direction x
T2 = A * x[0]   # Composante y du champ de torsion / Contrainte de cisaillement dans la direction y
T3 = ufl.as_ufl(0.0)  # Composante z nulle / Pas de contrainte dans la direction z (axe de torsion)

# Combinaison en un champ vectoriel de torsion
T = ufl.as_vector([T1, T2, T3])
print("Vecteur de contrainte de Cauchy défini pour la torsion (N/m²)")

Densité de contrainte de torsion (A_τ) = 2.7925e+09 N/m³
Vecteur de contrainte de Cauchy défini pour la torsion (N/m²)

# Définition des tenseurs de déformation et de contrainte
def epsilon(u):
    """Tenseur de déformation linéarisé"""
    return ufl.sym(ufl.grad(u))

def sigma(u):
    """Tenseur des contraintes (loi de Hooke)"""
    return lambda_ * ufl.nabla_div(u) * ufl.Identity(len(u)) + 2 * mu * epsilon(u)

# Définition des fonctions d'essai et de test
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)

# Définition de la force volumique (nulle dans ce cas)
f = fem.Constant(domain, default_scalar_type((0, 0, 0)))

# Formulation variationnelle du problème
a = ufl.inner(sigma(u), epsilon(v)) * ufl.dx  # Forme bilinéaire

# Créez un widget Dropdown
surface_selector = widgets.Dropdown(
    options=[('Surface latérale Sl', 'Sl'), ('Surface inférieure S0', 'S0'), ('Surface supérieure Sh', 'Sh')],
    value='Sh',
    description='Surface:',
    disabled=False,
)

# Créez un widget de sortie pour afficher les messages
output = widgets.Output()

# Fonction pour mettre à jour la forme linéaire
def update_linear_form(change):
    global L
    selected_surface = change['new']
    # Mise à jour de la forme linéaire
    L = ufl.dot(f, v) * ufl.dx + ufl.dot(T, v) * ds(surface_ids[selected_surface])
    # Mise à jour de l'affichage
    with output:
        output.clear_output(wait=True)
        print(f"Torsion appliquée sur {selected_surface}")

# Attacher la fonction de mise à jour au widget
surface_selector.observe(update_linear_form, names='value')

# Afficher le widget et la sortie
display(widgets.VBox([surface_selector, output]))

# Initialiser l'affichage
update_linear_form({'new': surface_selector.value})

# Créez un widget Dropdown pour la sélection de la condition aux limites
bc_selector = widgets.Dropdown(
    options=[
        ('Aucune', 'none'),
        ('Face inférieure S0', 'S0'),
        ('Face supérieure Sh', 'Sh'),
        ('Surface latérale Sl', 'Sl')
    ],
    value='S0',  # Valeur par défaut
    description='Encastrement:',
    disabled=False,
)

# Dictionnaire pour mapper les sélections aux conditions aux limites
bc_dict = {
    'none': [],
    'S0': [bc0_S0],
    'Sh': [bc0_Sh],
    'Sl': [bc0_Sl]
}

# Créez un widget de sortie pour afficher les messages
output = widgets.Output()

# Fonction pour mettre à jour le problème
def update_problem(change):
    global problem
    selected_bc = change['new']
    problem = LinearProblem(a, L, bcs=bc_dict[selected_bc], 
                            petsc_options={"ksp_type": "cg", "pc_type": "jacobi"})
    
    # Mise à jour de l'affichage
    with output:
        output.clear_output(wait=True)
        print(f"Encastrement appliquée sur {selected_bc}")
            
# Attacher la fonction de mise à jour au widget
bc_selector.observe(update_problem, names='value')

# Afficher le widget
display(bc_selector, output)

# Initialisation du problème avec la valeur par défaut
update_problem({'new': bc_selector.value})

uh = problem.solve()

# Création du maillage pour PyVista basé sur les coordonnées des dofs
u_topology, u_cell_types, u_geometry = plot.vtk_mesh(domain)

# Créez la grille PyVista et ajoutez les valeurs des dofs à la grille
u_grid = pv.UnstructuredGrid(u_topology, u_cell_types, u_geometry)

# Attach vector values to grid and warp grid by vector
u_grid["u"] = uh.x.array.reshape((u_geometry.shape[0], 3))

# Visualisation
p = pv.Plotter()
p.add_mesh(u_grid, show_edges=False, scalar_bar_args={
    "title": "u",
    "title_font_size": 24,
    "label_font_size": 22,
    "shadow": False,
    "italic": True,
    "font_family": "arial",
    "vertical": False
})
p.add_text("Déplacements", font_size=12, color="black", position="upper_edge")
p.add_axes(color="black")
p.set_background("grey")
p.show()

# Utiliser Panel pour créer un rendu interactif
pn.extension('vtk')
interactive_panel = pn.pane.VTK(p.ren_win, width=600, height=400)
from myst_nb import glue
glue("img", interactive_panel)

Visualisation des rotations principales

def omega(u):
    """Calcule le tenseur de rotation (partie antisymétrique du gradient de déplacement)"""
    return 0.5 * (ufl.grad(u).T - ufl.grad(u))

# Calcul du tenseur de rotation
theta = omega(uh)

# Extraction des composantes du tenseur de rotation
theta_xy, theta_xz, theta_yz = theta[0, 1], theta[0, 2], theta[1, 2]

# Définition de l'espace fonctionnel pour les rotations (éléments discontinus d'ordre 0)
V_omega = fem.functionspace(domain, ("DG", 0))

def interpolate_rotation(rotation_component):
    """Interpole une composante de rotation sur l'espace fonctionnel"""
    expr = fem.Expression(rotation_component, V_omega.element.interpolation_points())
    rotation = fem.Function(V_omega)
    rotation.interpolate(expr)
    return rotation

# Interpolation des composantes de rotation
rotation_x = interpolate_rotation(theta_yz)
rotation_y = interpolate_rotation(-theta_xz)  # Note le signe négatif
rotation_z = interpolate_rotation(theta_xy)

# Calcul de la norme de rotation
norm_theta = ufl.sqrt(theta_xy**2 + theta_xz**2 + theta_yz**2)
rotation_norm = interpolate_rotation(norm_theta)

# Configuration de la visualisation
pl = pv.Plotter(shape=(2, 2))
dargs = dict(cmap="jet", show_scalar_bar=False)
warped = u_grid.warp_by_vector("u", factor=0)  # Pas de déformation pour visualiser les déformations

def add_subplot(row, col, data, title):
    """Ajoute un sous-plot pour une composante de rotation"""
    pl.subplot(row, col)
    warped.cell_data[title] = data.vector.array
    warped.set_active_scalars(title)
    pl.add_mesh(warped, **dargs)
    pl.add_axes(color="k")
    pl.add_text(title, color='k', font_size=12)

# Création des sous-plots
add_subplot(0, 0, rotation_norm, "Norme de rotation")
add_subplot(0, 1, rotation_x, "Rotation RX")
add_subplot(1, 0, rotation_y, "Rotation RY")
add_subplot(1, 1, rotation_z, "Rotation RZ")

# Configuration finale
pl.link_views()
pl.camera_position = 'iso'
pl.background_color = 'grey'

# Ajout d'une barre de couleur globale
pl.add_scalar_bar(title="Rotation (rad)", n_labels=5, label_font_size=10, title_font_size=12, 
                  position_x=0.25, position_y=0.05, width=0.5)

# Affichage
pl.show()

Visualisation des déplacements amplifiés

# Création d'un objet Plotter PyVista pour la visualisation 3D
p = pv.Plotter()

# Application de la déformation au maillage
# Le facteur 0.04 amplifie la déformation pour une meilleure visibilité
warped = u_grid.warp_by_vector("u", factor=0.07)

# Ajout du maillage déformé à la scène
p.add_mesh(warped, 
           show_edges=False,  # Affiche les arêtes du maillage
           scalar_bar_args={
               "title": "Déplacement",  # Titre de la barre de couleur
               "title_font_size": 24,
               "label_font_size": 22,
               #"shadow": True,
               "italic": True,
               "font_family": "arial",
               "vertical": False,  # Orientation horizontale de la barre de couleur
               "n_labels": 5,  # Nombre d'étiquettes sur la barre de couleur
               #"fmt": "%.2e"  # Format scientifique pour les valeurs
           })

# Ajout d'un titre à la visualisation
p.add_text("Déplacements amplifiés", font_size=12, color="k", position="upper_edge")

# Ajout des axes de coordonnées
p.add_axes(xlabel='X', ylabel='Y', zlabel='Z', color="k")

# Configuration de l'arrière-plan et de l'éclairage
p.set_background("grey")
#p.enable_shadows()  # Ajoute des ombres pour améliorer la perception 3D

# Configuration de la caméra pour une vue optimale
p.camera_position = [(3, 3, 2), (0, 0, 0.5), (0, 0, 1)]

# Affichage de la visualisation
p.show()

Visualisation des déplacements normalisés et selon les axes principaux

# Définition des arguments communs pour l'affichage des maillages
dargs = dict(
    scalars="u",       # Utilise le champ de déplacement 'u' pour la coloration
    cmap="jet",        # Utilise la palette de couleurs 'jet'
    show_scalar_bar=False,  # Ne pas afficher la barre de couleur pour chaque sous-plot
)

# Création d'un plotter PyVista avec une grille 2x2 de sous-plots
pl = pv.Plotter(shape=(2, 2))

# Sous-plot 0,0 : Déplacement normalisé (magnitude)
pl.subplot(0, 0)
pl.add_mesh(u_grid, **dargs)
pl.add_axes(color="black")
pl.add_text("Normalized Displacement", color='k', font_size=10)

# Sous-plot 0,1 : Déplacement selon X
pl.subplot(0, 1)
pl.add_mesh(u_grid.copy(), component=0, **dargs)
pl.add_axes(color="black")
pl.add_text("X Displacement", color='k', font_size=10)

# Sous-plot 1,0 : Déplacement selon Y
pl.subplot(1, 0)
pl.add_mesh(u_grid.copy(), component=1, **dargs)
pl.add_axes(color="black")
pl.add_text("Y Displacement", color='k', font_size=10)

# Sous-plot 1,1 : Déplacement selon Z
pl.subplot(1, 1)
pl.add_mesh(u_grid.copy(), component=2, **dargs)
pl.add_axes(color="black")
pl.add_text("Z Displacement", color='k', font_size=10)

# Configuration globale de la visualisation
pl.link_views()  # Lie les vues des sous-plots pour une navigation synchronisée
pl.camera_position = 'iso'  # Vue isométrique
pl.background_color = 'grey'  # Couleur de fond grise

# Ajout d'une barre de couleur globale
pl.add_scalar_bar(title="Displacement (m)", n_labels=5, label_font_size=10, title_font_size=12, position_x=0.25, position_y=0.05, width=0.5)

# Affichage de la visualisation
pl.show()

Visualisation des composantes du tenseur des déformations

# Calcul du tenseur de déformation
eps = epsilon(uh)

# Extraction des composantes du tenseur de déformation
eps_xx, eps_xy, eps_xz = eps[0, 0], eps[0, 1], eps[0, 2]
eps_yy, eps_yz = eps[1, 1], eps[1, 2]
eps_zz = eps[2, 2]

# Création de la grille PyVista
eps_topology, eps_cell_types, eps_geometry = dolfinx.plot.vtk_mesh(domain)
eps_grid = pv.UnstructuredGrid(eps_topology, eps_cell_types, eps_geometry)

# Définition de l'espace fonctionnel pour les déformations (éléments discontinus d'ordre 0)
V_eps = fem.functionspace(domain, ("DG", 0))

# Fonction pour interpoler les composantes de déformation
def interpolate_strain(eps_component):
    expr = fem.Expression(eps_component, V_eps.element.interpolation_points())
    strain = fem.Function(V_eps)
    strain.interpolate(expr)
    return strain

# Interpolation des composantes de déformation
strain_components = {
    "XX": interpolate_strain(eps_xx),
    "YY": interpolate_strain(eps_yy),
    "ZZ": interpolate_strain(eps_zz),
    "YZ": interpolate_strain(eps_yz),
    "XZ": interpolate_strain(eps_xz),
    "XY": interpolate_strain(eps_xy)
}

# Configuration de la visualisation
pl = pv.Plotter(shape=(2, 3))
warped = u_grid.warp_by_vector("u", factor=0)  # Pas de déformation pour visualiser les déformations

dargs = dict(cmap="jet", show_scalar_bar=False)

# Fonction pour ajouter un sous-plot
def add_subplot(row, col, component):
    pl.subplot(row, col)
    warped.cell_data[f"Epsilon_{component}"] = strain_components[component].vector.array
    warped.set_active_scalars(f"Epsilon_{component}")
    pl.add_mesh(warped, **dargs)    
    pl.add_axes(color="black")
    pl.add_text(f"Epsilon {component}", color='k', font_size=12)

# Création des sous-plots pour chaque composante
for i, component in enumerate(["XX", "YY", "ZZ", "YZ", "XZ", "XY"]):
    add_subplot(i // 3, i % 3, component)

# Configuration finale
pl.link_views()
pl.camera_position = 'iso'
pl.background_color = 'grey'

# Ajout d'une barre de couleur globale
pl.add_scalar_bar(title="Strain", n_labels=5, label_font_size=10, title_font_size=12, position_x=0.25, position_y=0.05, width=0.5)

# Affichage
pl.show()

Visualisation des composantes du tenseur des contraintes

# Calculer le tenseur des contraintes
sig = sigma(uh)  # sigma est une fonction définie précédemment pour calculer les contraintes

# Extraire les composantes du tenseur des contraintes
sig_xx, sig_xy, sig_xz = sig[0, 0], sig[0, 1], sig[0, 2]
sig_yy, sig_yz = sig[1, 1], sig[1, 2]
sig_zz = sig[2, 2]

# Définir l'espace fonctionnel pour les contraintes (éléments discontinus d'ordre 0)
V_sig = fem.functionspace(domain, ("DG", 0))

# Fonction pour interpoler les composantes de contrainte
def interpolate_stress(sig_component):
    expr = fem.Expression(sig_component, V_sig.element.interpolation_points())
    stress = fem.Function(V_sig)
    stress.interpolate(expr)
    return stress

# Interpoler chaque composante de contrainte
stress_xx = interpolate_stress(sig_xx)
stress_yy = interpolate_stress(sig_yy)
stress_zz = interpolate_stress(sig_zz)
stress_yz = interpolate_stress(sig_yz)
stress_xz = interpolate_stress(sig_xz)
stress_xy = interpolate_stress(sig_xy)

# Configuration de la visualisation
pl = pv.Plotter(shape=(2, 3))  # Créer un plotter avec une grille 2x3

dargs = dict(
    cmap="jet",
    show_scalar_bar=False,
)

# Fonction pour ajouter un sous-plot
def add_subplot(row, col, stress_component, title):
    pl.subplot(row, col)
    warped.cell_data[title] = stress_component.vector.array
    warped.set_active_scalars(title)
    pl.add_mesh(warped, **dargs)    
    pl.add_axes(color="k")
    pl.add_text(title, color='k', font_size=12)

# Ajouter chaque composante de contrainte comme un sous-plot
add_subplot(0, 0, stress_xx, "Sigma XX")
add_subplot(0, 1, stress_yy, "Sigma YY")
add_subplot(0, 2, stress_zz, "Sigma ZZ")
add_subplot(1, 0, stress_yz, "Sigma YZ")
add_subplot(1, 1, stress_xz, "Sigma XZ")
add_subplot(1, 2, stress_xy, "Sigma XY")

# Configuration finale de la visualisation
pl.link_views()  # Lier les vues pour une navigation synchronisée
pl.camera_position = 'iso'  # Vue isométrique
pl.background_color = 'grey'  # Couleur de fond grise

# Ajout d'une barre de couleur globale
pl.add_scalar_bar(title="Contrainte (Pa)", n_labels=5, label_font_size=10, title_font_size=12, position_x=0.25, position_y=0.05, width=0.5)

# Afficher la visualisation
pl.show()

Calcul des contraintes de von Mises

s = sigma(uh) - 1. / 3 * ufl.tr(sigma(uh)) * ufl.Identity(len(uh))
von_Mises = ufl.sqrt(3. / 2 * ufl.inner(s, s))

V_von_mises = fem.functionspace(domain, ("DG", 0))
stress_expr = fem.Expression(von_Mises, V_von_mises.element.interpolation_points())
stresses = fem.Function(V_von_mises)
stresses.interpolate(stress_expr)

Visualisation des contraintes de von Mises

# Créez une nouvelle instance de Plotter
# Ajout des données de contrainte de von Mises au maillage déformé
warped.cell_data["VonMises"] = stresses.vector.array
# Définition de "VonMises" comme scalaire actif pour la visualisation
warped.set_active_scalars("VonMises")

# Création d'une nouvelle instance de Plotter PyVista
p = pv.Plotter()

# Ajout du maillage déformé avec la contrainte de von Mises
p.add_mesh(warped, 
           cmap="jet",  # Utilisation de la carte de couleurs "jet"
           show_edges=False,  # Ne pas afficher les arêtes du maillage
           scalar_bar_args={
               "title": "von Mises (MPa)",  # Titre de la barre de couleur
               "title_font_size": 24,  # Taille de la police du titre
               "label_font_size": 22,  # Taille de la police des étiquettes
               "shadow": True,  # Ajout d'une ombre à la barre de couleur
               "italic": True,  # Texte en italique
               "font_family": "arial",  # Police Arial
               "vertical": False  # Barre de couleur horizontale
           })

# Ajout d'un titre à la visualisation
p.add_text("Contraintes de von Mises", font_size=12, color="black", position="upper_edge")

# Ajout des axes de coordonnées en noir
p.add_axes(color="black")

# Configuration de l'arrière-plan en gris
p.set_background("grey")

# Configuration de la caméra pour une vue optimale
p.camera_position = [(3, 3, 2), (0, 0, 0.5), (0, 0, 1)]

# Affichage de la visualisation
p.show()

# Utiliser Panel pour créer un rendu interactif
pn.extension('vtk')
interactive_panel = pn.pane.VTK(p.ren_win, width=600, height=400)
from myst_nb import glue
glue("img", interactive_panel)

Analyses des résultats

Calcul des valeurs maximales

# Calcul du moment d'inertie polaire J pour une section circulaire
J = pi * (R ** 4) / 2

# Calcul de l'angle de torsion théorique
theta_theoretical = C * h / (mu * J)

# Extraction de la valeur maximale de la rotation RZ
max_rotation_z = np.max(np.abs(rotation_z.x.array))

# Extraction de la valeur maximale de la norme de rotation
max_rotation_norm = np.max(rotation_norm.x.array)

# Conversion en degrés pour une interprétation plus intuitive
max_rotation_z_degrees = np.degrees(max_rotation_z)

# Calcul de l'angle de torsion théorique
theoretical_twist = ang_radians * h  # ang est l'angle de torsion par unité de longueur, h est la hauteur du cylindre

# Comparaison pour RZ
relative_error_Z = abs(max_rotation_z - theoretical_twist) / theoretical_twist * 100

# Comparaison pour norme R
relative_error_N = (max_rotation_norm - theta_theoretical) / theta_theoretical * 100

print(np.max(rotation_z.x.array))

0.03471230845573148

# Calcul du moment d'inertie polaire J pour une section circulaire
J = pi * (R ** 4) / 2

# Calcul de l'angle de torsion théorique
theta_theoretical = C * h / (mu * J)

# Affichage des résultats
print(f"Moment d'inertie polaire J : {J:.6e} m^4")
print(f"Angle de torsion théorique θ : {theta_theoretical:.6f} radians")
print(f"Angle de torsion théorique θ : {np.degrees(theta_theoretical):.2f} degrés")

# Comparaison avec la rotation maximale calculée numériquement
if 'max_rotation_z' in locals():
    relative_error = abs(max_rotation_z - theta_theoretical) / theta_theoretical * 100
    print(f"\nComparaison avec la simulation numérique en RZ:")
    print(f"Rotation Z maximale simulée : {max_rotation_z:.6f} radians")
    print(f"Rotation Z maximale simulée : {np.degrees(max_rotation_z):.2f} degrés")
    print(f"Erreur relative : {relative_error_Z:.2f}%")
    
    print(f"\nComparaison avec la simulation numérique en norme:")
    print(f"Norme de rotation maximale : {max_rotation_norm:.6f} radians")
    print(f"Norme de rotation maximale : {np.degrees(max_rotation_norm):.2f} degrés")
    print(f"Erreur relative : {relative_error_N:.2f}%")
else:
    print("\nAttention : La rotation maximale simulée n'a pas été calculée précédemment.")

# Analyse des autres composantes
max_rotation_x = np.max(np.abs(rotation_x.x.array))
max_rotation_y = np.max(np.abs(rotation_y.x.array))

print("\nValeurs maximales des composantes de rotation :")
print(f"RX max : {np.degrees(max_rotation_x):.2f} degrés")
print(f"RX max  : {max_rotation_x:.6f} radians")
print(f"RY max : {np.degrees(max_rotation_y):.2f} degrés")
print(f"RY max  : {max_rotation_y:.6f} radians")

Moment d'inertie polaire J : 2.513274e-03 m^4
Angle de torsion théorique θ : 0.034907 radians
Angle de torsion théorique θ : 2.00 degrés

Comparaison avec la simulation numérique en RZ:
Rotation Z maximale simulée : 0.034712 radians
Rotation Z maximale simulée : 1.99 degrés
Erreur relative : 0.56%

Comparaison avec la simulation numérique en norme:
Norme de rotation maximale : 0.034734 radians
Norme de rotation maximale : 1.99 degrés
Erreur relative : -0.49%

Valeurs maximales des composantes de rotation :
RX max : 0.22 degrés
RX max  : 0.003798 radians
RY max : 0.21 degrés
RY max  : 0.003721 radians

# Extraction des valeurs maximales de rotation pour chaque composante
max_rotation_x = np.max(np.abs(rotation_x.x.array))
max_rotation_y = np.max(np.abs(rotation_y.x.array))
max_rotation_z = np.max(np.abs(rotation_z.x.array))

# Affichage des résultats
print("Valeurs maximales des composantes de rotation :")
print(f"RX max : {max_rotation_x:.6f} radians ({np.degrees(max_rotation_x):.2f} degrés)")
print(f"RY max : {max_rotation_y:.6f} radians ({np.degrees(max_rotation_y):.2f} degrés)")
print(f"RZ max : {max_rotation_z:.6f} radians ({np.degrees(max_rotation_z):.2f} degrés)")

# Calcul et affichage de la norme de rotation maximale
max_rotation_norm = np.sqrt(max_rotation_x**2 + max_rotation_y**2 + max_rotation_z**2)
print(f"\nNorme de rotation maximale : {max_rotation_norm:.6f} radians ({np.degrees(max_rotation_norm):.2f} degrés)")

# Comparaison avec l'angle de torsion théorique (si disponible)
if 'theta_theoretical' in locals():
    print(f"\nAngle de torsion théorique : {theta_theoretical:.6f} radians ({np.degrees(theta_theoretical):.2f} degrés)")
    
    # Calcul de l'erreur relative
    relative_error = (max_rotation_z - theta_theoretical) / theta_theoretical * 100
    print(f"Erreur relative (RZ max vs théorie) : {relative_error:.2f}%")

# Analyse de la contribution de chaque composante
total_rotation = max_rotation_x + max_rotation_y + max_rotation_z
print("\nContribution relative de chaque composante à la rotation totale :")
print(f"RX : {max_rotation_x / total_rotation * 100:.2f}%")
print(f"RY : {max_rotation_y / total_rotation * 100:.2f}%")
print(f"RZ : {max_rotation_z / total_rotation * 100:.2f}%")

Valeurs maximales des composantes de rotation :
RX max : 0.003798 radians (0.22 degrés)
RY max : 0.003721 radians (0.21 degrés)
RZ max : 0.034712 radians (1.99 degrés)

Norme de rotation maximale : 0.035117 radians (2.01 degrés)

Angle de torsion théorique : 0.034907 radians (2.00 degrés)
Erreur relative (RZ max vs théorie) : -0.56%

Contribution relative de chaque composante à la rotation totale :
RX : 8.99%
RY : 8.81%
RZ : 82.20%

import numpy as np

# Valeur analytique de la rotation maximale (ici, 0.035 radians)
theta_analytique = 0.035

# Valeur numérique récupérée
theta_numerique = max_rotation_z

# Calcul de l'erreur absolue
erreur_absolue = np.abs(theta_analytique - theta_numerique)

# Calcul de l'erreur relative
erreur_relative = erreur_absolue / np.abs(theta_analytique) * 100

# Affichage des résultats
print(f"Rotation analytique maximale : {theta_analytique} radians")
print(f"Rotation numérique maximale : {theta_numerique} radians")
print(f"Erreur absolue : {erreur_absolue} radians")
print(f"Erreur relative : {erreur_relative} %")

Rotation analytique maximale : 0.035 radians
Rotation numérique maximale : 0.03471230845573148 radians
Erreur absolue : 0.0002876915442685257 radians
Erreur relative : 0.8219758407672162 %